Stochastic unraveling of Redfield master equations 
and its application to electron transfer problems 

Ivan Kondov/'Q Ulrich Kleinekathofer, 2 and Michael Schreiber 1 

1 Institut fur Physik, Technische Universitat, 09107 Chemnitz, Germany 
2 International University Bremen, P.O.Box 750 561, 28725 Bremen, Germany 

Abstract 

A method for stochastic unraveling of general time-local quantum master equations (QMEs) is 
proposed. The present kind of jump algorithm allows a numerically efficient treatment of QMEs 
which are not in Lindblad form, i.e. are not positive semidefinite by definition. The unraveling can 
be achieved by allowing for trajectories with negative weights. Such a property is necessary, e.g. 
to unravel the Redfield QME and to treat various related problems with high numerical efficiency. 
The method is successfully tested on the damped harmonic oscillator and on electron transfer 
models including one and two reaction coordinates. The obtained results are compared to those 
from a direct propagation of the reduced density matrix (RDM) as well as from the standard 
quantum jump method. Comparison of the numerical efficiency is performed considering both the 
population dynamics and the RDM in the Wigner phase space representation. 
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I. INTRODUCTION 



Time-independent as well as time-dependent phenomena in chemical physics, quantum 
optics, solid state physics, biological physics, etc. are often described using QMEs 1 ? 2 . In 
particular, electron transfer (ET) dynamics in open quantum systems, i.e., systems with 
dissipation, can be conveniently treated within this formalism. The QMEs govern the time 
evolution of density matrices which are used in order to represent the mixed nature of the 
states. Recently stochastic wave function methods have received a great deal of attention. 
In unraveling schemes one considers an ensemble of stochastic Schrodinger equations (SSEs) 
which in the limit of a large ensemble resemble the respective QME. Although all SSE 
approaches have a common basis 3 they are usually divided into two classes. One is the 
quantum jump method also known as Monte Carlo wave function (MCWF) approac h 4 ' 5 ' 6 '?! 8 . 
In this approach the dynamics is described by a Schrodinger-like wave equation interrupted 
by instantaneous deviations from the continuous motion (quantum jumps). The second class 
of SSE approaches are the quantum diffusion models with continuous motion^ifi which are 
not in the center of interest here. The numerical effort for solving SSEs scales much more 
favorably with the size of the basis than a direct propagation of a density matrix since one is 
dealing with wave functions and not density matrices (for a comparison of direct integrators, 
see Ref. lilt). Thus, stochastic unraveling is an efficient numerical tool for solving QMEs. Of 
course, to achieve good statistics one has to average over a large number of wave functions. 
So the SSE approaches become preferable for large and complex systems with many degrees 
of freedom. In passing, we want to mention that in the present paper we are not interested 
in establishing a relation between the SSE dynamics and some measurement process. Thus, 
the treatment of single trajectories will be done without giving a special physical meaning 
to them. 

One of the important properties of a density matrix is the positive semi-definiteness for 
all times, i.e. that all populations are positive or zero. This property is fulfilled for QMEs 
of the Lindblad formi 2 - but not necessarily for reduced dynamics in genera l 13 ' 14 . A slightly 
generalized generator for a completely positive density-matrix evolution can be found in 
Ref. 15 while a discussion on the non-Markovian case has been done in Ref. 0. Most of the 
unraveling schemes^ ' 6 ' 9 ' 10 ' 17 ' 18 have been restricted to QMEs of Lindblad formi 2 - which ensures 
that the RDM stays positive semidefinite for all times and all parameters. Nevertheless, 
there are many physically meaningful QMEs which result in positive semidefinite RDMs 
although they are not of Lindblad for m 19 ' 2 *- 1 . The increasing interest in descriptions beyond 
the Lindblad class such as the quantum Brownian motio n 19 ' 21 , the Redfield formalism 2 * 2 ^, 
non-Markovian scheme o 2 ^' 24 ' 25 , etc. resulted in various efforts to develop new stochastic wave 
function algorithms. 

Strunz et a/ . 19 ' 21 extended the QME for Brownian motion to a non-Markovian QME and 
then applied a quantum state diffusion algorithm. A similar approach was also proposed 
by Gaspard et aZ.— . Recently Stockburger and Graberl^ developed a method for an exact 
formulation of the RDM in terms of SSEs of a system coupled to a linear heat bath. Breuer et 
ali& extended a scheme which they had used to calculate multi-time correlation functions^ 
to the unraveling of QMEs. Their technique is based on doubling the Hilbert space. So 
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instead of a single stochastic wave function one has a pair of them 25 . In this approach, 
norm and Hermiticity are not preserved in single realizations but only in the ensemble 
average which makes the algorithm unstable. Since stability and efficiency are crucial issues 
for unraveling algorithms we propose an alternative approach which fulfills these criteria. 
Though the present approach has only been tested for Redfield and Brownian-dynamics 
master equations so far 2 - there are only few restrictions to its range of validity and it is 
therefore applicable to a much larger class of time-local quantum master equations. 

ET is commonly treated in modern theories with use of the RDM formalism and QMEs. 
Alternatives to QMEs are, for example, semi-classical theories^, path integral method o 31 ' 32 
and recently the self-consistent hybrid approac h 33 ' 34 ' 35 ' . The latter was shown to treat suc- 
cessfully the spin-boson problem 3 ^, the ET in mixed-valence compound o 34 ' 35 as well as the 
heterogeneous ET at semiconductor surfaces 35 . Solving the QME for ET in model systems 
with on o 36 i 37 i 38 i 39 and many reaction coordinates 4 ^ has been done with success. Exhaustive 



reviews on ET can be found, e.g., in Refs. |4l| and|42l Apart from the non-Markovian de- 
scriptions of transfer phenomen a 23 ' 43 ' 44 ' 45 the use of Redfield theory for ET was investigated 
as wel l 11 ' 37 ' 38 ' 46 ' 47 ' 48 ' 49 . The model used in the latter references is based on vibronically cou- 
pled diabatic potentials which are sufficiently well approximated by harmonic potentials. In 
particular, the influence of the electronic coupling between the diabatic states on the dissipa- 
tion was investigated. Neglecting this effect results in the diabatic damping approximation 
(DDA) 4 ^ 47 . This approximation as well as considering first order perturbation theory in 
the electronic coupling were objects of recent studies^ 7 -. A typical problem that occurs with 
increasing the complexity of the ET models, i.e. the dimension of the RDM, is the numerical 
effort. Thus, the stochastic unraveling of generalized time-local QME was developed with 
the prospect of applications to more complex ET systems. 

Recently the present scheme was briefly demonstrated for the quantum Brownian motion 
of a harmonic oscillator 29 . In the present paper the stochastic unraveling of the Redfield 
QME shall be considered in more detail as well as applications concerning multi-mode models 
for ET shall be presented. In the next Section a brief introduction to the Redfield formalism 
will be given. Section ITTT1 focuses on the derivation of the SSEs relevant for the generalized 
time-local QME while Section IIVI will provide explicit expressions for the jump rates. In 
Section |V| we describe three concrete applications of the proposed quantum jump method: 
the damped harmonic oscillator and a model for ET with one and two reaction modes. 
A study and discussion of the numerical efficiency in Section IVII and a conclusion follow. 
The detailed quantum jump algorithm used in the present contribution can be found in the 
Appendix. Atomic units are used throughout the paper, i.e. h = 1. 

II. REDFIELD FORMALISM 

In Redfield theory the overall system is partitioned into a relevant system whose evolution 
is of interest and a thermal bath using the Hamiltonian 

H = H S + H B + H SB . (1) 
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Here H$ and H-q are the Hamiltonians of the relevant system and of the bath while i^sB 
describes their interaction. In Subsec. IV CI it will be shown how this partitioning can be 
rigorously performed. In general, the interaction part can be represented by bilinear products 

#sb = ® mKm ( 2 ) 

m 

of system and bath operators, K m and $ m , respectively. In the following K m and $ m will 
be considered Hermitian. The state of the system is described by the RDM performing a 
trace of the total density matrix a over the bath degrees of freedom, i.e. p = tv^a. It is 
assumed that the bath stays in thermodynamic equilibrium at all times. This means that 
the relaxation of the bath is much faster than the evolution of the system. In addition, one 
assumes that the system-bath interaction is sufficiently small to be treated perturbatively to 
second order. Using the Hamiltonian (0) and the assumptions described above one obtains 
a non-Markovian QME for the RDM. One possible way to obtain a Markovian QME instead 
is to neglect memory effects which are due to the finite bath correlation time. The formal 
treatment yields the Redfield QMB^i 

p=-i [H 8 , p] + { [ A ^' K ™ ] + l K ™ P A U } • ( 3 ) 

Note that the operator h^ m is the adjoint of the relaxation operator A m . This is only true if 
K m and $ m are Hermitian^ii but not in genera l 22 ^ . The relaxation operator is given by 

oo 

A m = J2 f drC mn {r)K l n {-r), (4) 

n 

where C mn (r) is the bath correlation function and K\ the system operator in the interaction 
picture. Usually it is easier to obtain the latter quantity in the frequency domain, e.g. with 
use of molecular dynamics simulations or with a simple bath modeling. Either approach 
yields the bath spectral density J mn (oj) in terms of which the correlation function can be 
constructed as 2 - 

C mn {u) = 2tt[1 + n(u})][J mn (u) - J mn {-U))] (5) 

with the Bose-Einstein distribution n(u) = (er* — and (3 = (ksT) -1 being the inverse 
temperature. All considerations in the present work will be limited to the Ohmic form of 
the spectral density with exponential cut-off. However, a spectral density of Debye form can 
be constructed which results in nearly the same values of J mn (uj) for the specific spectrum 
of H$ used. 

Under certain approximations Eq. Q can be transformed to Lindblad forrr4& 



dp(t) 
dt 



-i [H s ,p(t)] + J2 



LnP(t)Ll - l -p(t)LiL n - hiL n p(t) 



(6) 



One way to obtain the Lindblad QME (jUJ) is starting either from the non-Markovian QME 
or from Eq. (jlj) and assuming a 5-correlated bath 2 ^, i.e., C mri (r) — > c mn 6(r). Subsequent 



4 



diagonalization of the correlation matrix c by means of a unitary transformation c = U^kU, 
yields Eq. (jHJ) with L n = -^/K^^2 m U nm K m . Alternatively, the Lindblad QME is obtained 
when the DD A 46 ' 47 is invoked with the so-called rotating wave approximation (RWA) in the 
system-bath coupling as done in the present paper. The explicit form of the corresponding 
Lindblad operators used will be given in the sections below. 

The Lindblad QME (JHJ) will be used to compare the numerical efficiency of the present 
quantum jump method method with the standard on o 4 '^' 1 ^ 1 !' by solving the single-mode 



and the two-mode ET model. For more details of the DDA we refer to Refs. |4y and |42 



III. STOCHASTIC SCHRODINGER EQUATION 

To start with the derivation of the unraveling scheme we first like to formulate the time- 
local QME in its most general form 

^ = A(t)p(t)+p(t)A\t) 

AI 

+ £{C fe (t)p(t)4(t) + E k (t)p(t)Cl(t)} . (7) 

k=l 

Because A(t), Ck(t), and E^it) are arbitrary operators this equation conserves only Her- 
miticity. In order to conserve also the norm further restrictions have to be applied as shown 
below. All time arguments will be dropped henceforth for clarity. 

The RDM will be recovered by averaging over an ensemble of two vectors \ip) and 
which are elements of the doubled Hilbert space, as 



P= 1^(01 + 10)^1 • (8) 

Every individual realization of the stochastic process before averaging denoted by the pair 
|0)) will be called a trajectory. In contrast to Ref.|25|the averaging formula (jSJ) preserves 
Hermiticity of single trajectories leading to a significantly improved numerical performance 
of the scheme. Each trajectory \<p)) is propagated be means of two SSEs having the 
following generic form 



M 2 



k=l i=l 
M 2 

d\<t>) = D 2 \<t>)dt+Y<Y< si ^ d & ■ ( 9b ) 



k=l i=l 



Unlike deterministic differential equations the SSEs include differentials of the complex noise 
variables in addition to the time variable. The superscript in Q. denotes which of the 
two terms from the Hermitian pair in the sum in Eq. (0) is taken while the subscript 
counts the relevant dissipative channels from 1 to M. The operators Di and D 2 specify the 
deterministic and the operators Sj k the stochastic part of the evolution. In general, they 
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may be time- dependent. The stochastic differentials^ d^ l k are assumed to have zero mean, 
to be uncorrelated and normalized to dt: 



da = 0, d&d£i = 5 tJ S kl dt . (10) 

Differentiating Eq. (jHJ), neglecting all terms higher than first order in dt, and assuming that 
ensemble averages always factorize^ yields 



dp = [^1^(01 + ^10)^1 



dt 



M 

+ E K WM s ll + s lk \WMs!l] dt + n.c. (ii) 



k=l 

Comparing Eq. with the original QME (JJJ) one is able to replace 

Slk = s lk = C k + a\ and S\ k = Sl k = E k + a\ (12) 

where a\ and ot\ are arbitrary possibly time-dependent scalar functions of (\ip), |0)). Plug- 
ging the latter expressions into Eq. (jTTj) yields 

M 

D 1 = D 2 = A-J2 H*C k + al*E k + ala 2 k *) . (13) 

k=l 

According to Ref. Eq. Q describes a quantum diffusion process if the leading terms in 
dC,l are of first order in y/di. When dQ, can be given by a finite number of values only, e.g. 
±\/di, the process results in continuous but random trajectories within each infinitesimal 
time interval dt (for dt — > the trajectories become smooth but still stay noisy). In that way 
one derives diffusion methods which will not be considered in the present work. However, if 
the leading terms in d£ l k have finite values of order unity, i.e. zeroth order in \/dt, Eq. leads 
to the so-called quantum jump methods which produce trajectories that are deterministic 
during finite time intervals connected by discontinuous transitions (jumps). The jumps are 
specified by their jump rates p k , which have to be real scalar functions of (|-0), |0)). If n\{t) 
is the number of jumps in channel k due to term i up to time t, the probability for n k (t) to 
increase by one, i.e. the expectation value of both dn l k and (dn k ) 2 , should be equal to p k dt 
during the infinitesimal time interval dt. This can be written as^ 

da = dnlk ~ft d * e * (14) 



so that it obeys condition (JTUJ). The phase factor e lLp leads merely to a phase shift in the 
wave vectors and cancels within each realization and we therefore set ip = 0. If dn k vanishes 
for all k and i, then Eq. Q becomes a deterministic Schrodinger equation. For any k 
and i, dn\ = 1 indicates the occurrence of a jump. In this case we have d(\<f>), = 
(10), m after jump (10), l^))beforejum P - Taking this into account and substituting Eqs. (JHJ) 
and (|12|) into Eq. Q it is found that a\ = —\fp\- Eventually, the final form of the SSEs 
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for the quantum jump method is obtained as 

M 



d\<j>) 



Pk + Pk 



k=l 



M 



El 



k=l 



M 



k=l 



\ip)dt 



— 1 dnl 



k=l 

M 1,2 



' -l\dnl 



— 1 )dn 2 k 



E, 



1 dnl 



(15a) 



(15b) 



IV. JUMP RATES 



Essential for the performance and particularly for the convergence behavior of the quan- 
tum jump method is how the jump rates p\ and p\ are specified. They have no physical 
meaning since in the average they do not influence any observable but determine the statis- 
tical error. A detailed discussion on the optimization of the jump rates as free parameters 
can be found in Ref. |52|. Another freedom is that Eq. (J7J) is invariant with respect to a 
gauge transformation of the kind Ck fCk, E^ — > Ej./ f if / is a real, scalar function of 
time. Each single realization, and hence the stochastic process, is independent of this gauge 
transformation and using such transformation offers us no further advantages. We note that 
the jump rates in Ref. |25j do not fulfill the invariance under this gauge transformation. 

Following the approach in Ref. we require that the norm of every single trajectory is 
constant in time. Under such a condition expressed as 



d 

u< it 







(16) 



the p k are adapted at each moment of time. This approach yields a numerically stable and 
efficient algorithm. In contrast, numerical tests with jump rates adapted to other quantities 
such as (0|0), e ^ c - resulted in an unstable scheme. The operators that enter the 

QME (JZJ) are restricted by condition (fTB*j) yielding 



M 

A + At + E + c l E k) = °- 



(17) 



k=l 



Let us try to determine the jump rates from this condition. The total jump rate is obtained 
applying Eq. (fIB)l to the deterministic part of Eq. (fT5|): 



\A + A^\^) + (^\A + A^\ 

<0|V) + (V|0> 



(18) 
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All partial jump rates can be then successively found using Eqs. (|T7|) and (fTHI) : 



Pi 



Pi 



\4c k \ 



\E\C k 



<0i-0> + (m 



(19a) 
(19b) 



Here a problem occurs because the values of the p\ do not have to be positive for all 
trajectories at all times. But, since the jump rates p\ are arbitrary real functions we can 
choose them as the absolute values of the p\ 



Pk 



Pk 




(20a) 
(20b) 



An additional weight factor ±1 for the trajectories has to be introduced which changes 
its sign every time a jump is performed with p\ = —p\. It can also be implemented (as 
in the appendix) by allowing for negative norms of the trajectories. The change from p l k 
to p k gives rise to a small deviation of the norm from unity because in the regions where 
not all p k and p\ are identical, norm conservation is no longer guaranteed, i.e. the sum of 
the p\ differs from p. As long as the occurrence of a jump is a very rare event and the 
number of negative p l k is also very small the deviation from the initial norm is expected 
to be small. In all tests this deviation was far below 1% without effecting the numerical 
efficiency of the proposed algorithm. This is how the scheme tolerates trajectories with 
possibly negative weights which arise from the fact that the RDM with the QME ((7j) is not 
necessarily positive semidefinite. If the RDM stays positive semidefinite during its entire 
time evolution the negative weights are not needed, i.e. all trajectories can be normalized 
to unity and represent physically relevant pure states. A possible implementation of the 
present unraveling scheme is shown in the appendix. 

In the examples below the RDM can exhibit negative eigenvalues. This nonphysical 
situation could probably be improved by applying an initial slippage to the initial stated 
Another possibility to avoid non-positive semidefinite RDMs is to start with a derivation 
of different QMEs in the form (J7|) with time-dependent coefficients. It has been shown 
that non-Lindblad QMEs with time-dependent coefficients can preserve the positivity of the 
RD M 19 i 20 . Nevertheless, an unraveling scheme has to be able to follow also the nonphysical 
behavior of the QME because in the ensemble average the solution of the SSEs should 
completely coincide with the exact solution of the QME. 
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V. APPLICATIONS OF THE PRESENT UNRAVELING METHOD 



A. Damped harmonic oscillator 



One of the most simple toy models used for testing in dissipative quantum dynamics is 
the damped harmonic oscillator. Here it will be formulated within Redfield theory, i.e. one 
has to obtain the explicit form of the operators in Eq. (jSJ). The oscillator has mass M. 
and frequency u . If the thermal bath is modeled by quantum harmonic oscillators and the 
system operator K is the oscillator coordinate q = (a + al) / y/2M.u)Q the sum in Eq. (J3J) 
contains only one term in which 



A = r^^lp [(n(u ) + 1) a + n(u )a^} . (21) 

The damping rate T is related to the spectral density of the bath J (no) as T = 7iJ(uj ) / (A4uj ) . 
Therefore, the explicit form of J{oS) is not necessary since the oscillator and the bath in- 
terchange quanta only at the frequency Uq. Performing either the RW A 53 1 54 or the secular 
approximatio n 2 ^ 4 '^ Eq. (J5J) is transformed into a Lindblad QME @. For the sake of simplic- 
ity this will be shown here with the RWA for the harmonic oscillator but the generalization 
for the ET model solved within the DDA is straightforward. Inserting the expressions for A 
and K into Eq. (JHJ), denoting b\ = a and b 2 = a , and performing some calculus the QME 
obtains the form 

p = -i [h q , P ] + ^ Yl K a ( b ^ b l - l b fap - lp b h) > ( 22 ) 

i,j=i ^ ' 

where k is the correlation matrix 



n = ( x n W + \ V (23) 

In order to transform Eq. (|22|) into Lindblad form either k, has to be diagonalized imposing 
conditions for which the eigenvalues are positive^ or the RWA in the system-bath coupling^ 
has to be performed. It is easily seen that the determinant of n is — 4 and hence the former 
method fails for this QME. Performing the RWA implies that the off-diagonal elements of 
k, are set to zero. Then the Lindblad operators take the explicit form 

L L = vS) + l)Ta and L 2 = v / n(o;o)ra t . (24) 

This result is easily generalized for the ET models discussed in the next two subsections. L\ 
and L2 have a clear physical interpretation. L\ damps all occupied levels bringing their pop- 
ulations one level lower, while L2 has the opposite effect. In the thermodynamic equilibrium 
the jump rates for both operators are equal: the populations do not change. 

To find the operators involved in the generalized QME (J2J) (with M — 1 for the harmonic 
oscillator) one has to carry out the commutators in Eq. Q. Then one can easily identify 

C X = K, Ei = A, A = -iH s - KA. (25) 
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In contrast to the Lindblad operators (|24|) the action of C\ and E% on the wave function 
is more subtle. This is why it is difficult, and probably not possible, to assign a certain 
physical process to a single trajectory. 

In our stochastic simulation of the damped harmonic oscillator the temperature T = ujq/A 
and T = uo/10 are used. Figure [T] shows the population dynamics of the lowest four levels 
of the oscillator starting from the pure initial state P33 = 1. As seen, the convergence to 
the exact solution is very slow (10 4 trajectories are still not sufficient). On the other hand 
the test system is very small and the QME can be solved very fast using direct propagators. 
The true advantage of the method can be seen with larger systems, where it shows both a 
faster convergence and a good scaling. 



B. Electron transfer model with one reaction coordinate 



Let us consider a model for electron transfer with the Hamiltonian^ 

H S = H^+V = Y^Hi\i){i\+V (26) 

i 

where Hi are the Hamiltonians of two harmonic oscillators (i.e., i = 1,2) which describe the 
vibronic spectrum of two electronic states interacting via the electronic coupling V. If the 
system includes a single reaction coordinate q the vibronic Hamiltonians read 

H i = U l + u)i{a\oi + ^) + ^(o* + 4), (27) 

where a« and aj are the boson operators, Aj the dimensionless displacements of the harmonic 
potentials along the reaction coordinate, Ui the oscillator frequencies, and U{ the electronic 
excitation energies. A very useful parameter of the system which is related to the last 
term in Eq. (|2*7jl is the reorganization energy Aj = cjjA?/2. It is also proportional to the 
vibronic coupling cjjAj. Using the former expression one can define the potential minima 
as Uf = Ui — Aj. Configurations in which the potential minimum of the upper free-energy 
surface is lower in energy than the lower free-energy surface at that point are in the so- 
called normal region. If the opposite is true the configuration is in the Marcus inverted 
region. The electronic coupling i>i 2 between the model potential surfaces is independent of 
the coordinate. So the respective term in Eq. (J26j) obtains the explicit form 

V = Y,^0--8iiHf(hM;j,N)\iM)(jN\ . (28) 

i,j M,N 

The Franck-Condon factors f(i, M;j, N) are calculated using the eigenfunctions cpiM of the 
harmonic oscillators 



f(i,M;j,N) = (%M\jN) = J dqip m {q)ip jN {q) . (29) 

By analogy with the damped harmonic oscillator the system operator K is defined as the 
coordinate operator, i.e. 

K = q = Y^ (2uiMy l/2 (at + a*) . (30) 
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We consider a potential configuration in the normal region with no barrier between the 
two harmonic potentials which have equal curvature {uo\ = cu 2 = ^o), change of free energy 
— U[ ^ = —2u , reorganization energy Ai = 3lu , A2 = and inter-center coupling 
v 12 = ujq. The reason for this choice is the intention to compare the standard and the new 
quantum jump methods. As the Lindblad QME is obtainable only with RWA and DDA 
one has to study a parameter region where both QMEs generate almost the same dynamics. 
The bath is described by an Ohmic spectral density with cut-off frequency uj c = u at 
temperature k^T = ujq/A. The system-bath interaction is characterized by the damping 
rate Y = 7r^/(A4 exp(l)) = cuo/10 (see Ref. 1471 for details). Again, the operators necessary 
for the new quantum jump method are defined according to Eq. (|25jl . 

A Gaussian wave packet located at the donor state |1) having energy slightly above the 
crossing of the harmonic potentials was chosen as initial state. The numerical simulation with 
about 1000 trajectories provides sufficiently converged and accurate results. Figure El shows 
the relaxation of the ensemble-averaged donor population P x = (t/>|1)(1|0) + (<f)\l) . 
A widely discussed property of the Redfield equation is that it does not strictly conserve 
positivity of the RDM£4 Although Pi is always positive the tiny negative fraction in Fig. El 
is an evidence for the existence of single realizations with negative P±. In contrast, the 
simulation of the same system with the Lindblad operators ()24j) by means of the standard 
quantum jump method^*^*^ keeps all values of Pi well confined between and 1. 

Besides numerical efficiency, another advantage of the quantum trajectories is the better 
insight into the quantum mechanisms underlying the overall dynamics of the ensemble. 
Though it is impossible to give a direct physical interpretation of every single trajectory 
one can extract information from the ensemble statistics. As we can see in Fig. El the 
distribution of the individual expectation values of the population is skew and comprises 
several maxima. One can better visualize the wave packet dynamics in phase space using 
the Wigner representation of the RDM^ as done in Figs. H] and El The evolution of the 
expectation values of the momentum and the coordinate, which can be regarded as the 
center of mass of the wave packet, is described by a path in phase space as shown in Fig. El 

The wave packet starts off with zero momentum from the location of the excited state. 
One can distinguish two stages of the ET dynamics. In the first stage the wave packet splits 
into several parts and occupies the whole accessible phase volume, i.e., it spreads (see Figs.|U 
b, c andEJb, c). The motion of the principal part of the wave packet, which is seen as a sharp 
peak in Fig. El implies coherent transfer of population. The peak moves rapidly with time 
in an oscillatory fashion while its amplitude decays as decoherence processes advance. After 
that all parts of the wave packet coalesce to a single bell-shaped distribution. In the second 
stage the wave packet continues to propagate slowly in phase space while its maximum is 
approaching the equilibrium point at (p)t-*oo an d (q)t->oo- This slow motion is seen as a 
small drift to the right beginning from the central region of the spiral path (see Fig. EJ). It is 
due to dissipative transfer mechanisms 4 - and is small for barrierless potential configurations 
as in the present case. We note also that in the time between the third population revival 
and the 10th vibrational period (see Fig. |2J one can recognize a crossover between the two 
stages discussed above. Such a moment of time is shown in Figs. El lU and El 

Unlike the trajectories considered in Ref. within the quantum diffusion approach the 
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wave packet of one individual trajectory in our calculation spreads over the whole phase 
volume of the system (Fig. @]a). One reason for this discrepancy is the different value used 
for h. Generally, for sufficiently small H the trajectories become classical states and virtually 
shrink to points in phase space. However, the problem of their localization for small % is 
non-trivial because the system may be almost classical but with chaotic behavior. 



C. Electron transfer model with multiple reaction coordinates 

Modeling the system-bath separation one has to minimize the degrees of freedom in the 
relevant system and the system-bath coupling simultaneously. The small polaron transfor- 
mation has been used to effectively reduce the system-bath coupling for a two-level spin- 
boson system^. It also has been shown that this approach can be extended for multi-level 
systems^. Alternatively, it is possible that one can successively take strongly coupled de- 
grees of freedom from the bath and put them into the relevant-system part. This will make 
the effective system-bath coupling smaller and hence the application of the Redfield theory 
more reasonable. Multi-mode modeling of ET reactions, including systems in the inverted 
region, has been done in Refs. 00, and with similar argumentation. On the other 



hand, there is experimental evidence for the participation of multiple modes in the ET tran- 
sition in some systems, such as oxazine-1 in N,N-dimethylaniline^ and betaine-30 in various 
solvents^. Correspondingly, the relevant part of the total ET system can be modeled with 
a treatable small set of R reaction coordinates {qi\. For this purpose one may select a set 
of representative harmonic normal modes from the molecule and from its environment (e.g. 
the solvent or the crystal lattice). Since all normal modes are decoupled one can use the 
single-mode operators to calculate the matrix elements of the necessary operators in the 
diabatic basis \iM\ . . . Mr). The Hamiltonian of each diabatic electronic state reads 



R r -j A- 

Hi = Ui + J2 ^,1(4,1^,1 + 2") + K* + 4,i) 



(31) 



For R reaction modes Eq. (J2J) includes R summation terms linear in each coordinate q\. The 
matrix element of Ki reads 

(iMt . . . M R \Ki\jN! ...N R )= — 6a (SMt+i,*, + 1 + 6 Mi -im v 7 ^) II <^4 32 ) 

where A4 is the reduced mass of the relevant system. 

The multi-mode model can be easily reduced to an effective single-mode model by means 
of an orthogonal transformatio n^ 1 of the Hamiltonian H^°\ For two diabatic states with 
equal curvatures one can drop the electronic index. Denoting the relative displacement by 
A; the transformation has the form^i 

R 

B k = J2 Vkl(Xl ' ( 33 ) 
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with 



Vfa = m A '"' „ , Dl = J2 A '"' for MO (34) 



and 



Vr 



R 



0/ 



(W . (W = £(A,.,) 2 , "o-^E-?A?. (35) 
For 7^ the new mode frequencies are the roots of the equation 

E^T = 0. (36) 
i=i 1 K 

After this transformation the multi-mode Hamiltonian is cast into the form (P) where 



tt A 

H B = Y,^B[B k (38) 



H s = U + n BlB + -^(A) + 4) (37) 



R-l 



k=l 

H SB = J2^(BlB h + BlB ). (39) 

k=l k 

One can see that all R normal modes are transformed to a finite bath with R — 1 modes 
with frequencies Qk- A new effective mode with frequency Qq and displacement Ao is 
created which is bilinearly coupled to the new bath modes. From Eq. ()35j) it follows that 
A Q < Aiu t for positive coordinate displacements, i.e. the vibronic coupling of the new 
effective mode is smaller than the total vibronic coupling of the initial normal modes. Thus, 
a reduction from a multi-mode to a single-mode model extends the bath and enlarges the 
system-bath coupling. This reduction is unique. On the other hand, the addition of a bath 
mode to the relevant system is not unique. It depends on the choice of a certain bath mode. 
Nevertheless, it always reduces the system-bath coupling and enlarges the total vibronic 
coupling of the relevant system. 

In the following, a two-mode ET model will be considered. The frequencies of the two 
modes were chosen 0.07 and 0.18 eV, and the reorganization energies 0.33 and 0.82 eV, 
respectively. These frequencies correspond to internal molecular modes in real ET systems. 
The effective mode was calculated using Eq. ()35j) . The change of the free energy was taken 
to be U$ — U® — —0.2 eV and the electronic coupling between the diabatic electronic states 
v\2 = 0.1 eV. Again, a harmonic bath with an Ohmic spectral density was considered 
with damping rate T = 0.007 eV and temperature 295 K. 16 levels for each mode gave 
a good convergence. Initially, the lowest vibrational level of the excited state is populated 
(Fig.[7|). The coherent dynamics of the excited state of the effective mode model yields small 
population in the ground state. This is followed by an almost complete revival at about 110 
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fs. In contrast, the two- mode model exhibits an ultrafast coherent population transfer to 
the ground state. Recurrences of population in the excited state appear many times within 
120 fs but their yield does not exceed 50 %. As the two-mode system can be regarded as a 
conservative two-particle system its complete recurrence period (Poincare's cycle) is much 
longer than that for the single-mode system. In a way the isolated two-mode system shows 
a relaxation behavior that is typical for open systems. 

Turning on the dissipation leads to irreversible transfer to the lower state. However, 
the picture does not change qualitatively as the dissipative transfer mechanism does not 
contribute significantly at the early times shown here. As already founo™^! and seen in 
Fig. the DDA has also no serious influence on the early dynamics. As the DDA has been 
recently well studied we have put the focus here on the ultrafast initial stage where the 
difference between the single-mode model and the two-mode model is best characterized. 

VI. EFFICIENCY AND STABILITY 

To estimate the convergence and hence the stability of the proposed scheme one needs 
an appropriate measure for the convergence. Unraveling the QME one aims to calculate 
the time evolution of an observable, e.g. the population P 1 , performing an average over N s 
single trajectories. At time tj the average reads 



where N t is the total number of propagation time steps. In this convergence measure we use 
as reference the average performed over N s — k trajectories. For convenience the increment 
k can be chosen to be the number of computing nodes in a parallel implementation of the 
stochastic algorithm. One can easily see that e vanishes for large N s if both terms in the 
sum converge to Pi(tj, N s — > oo). If these two terms diverge with increasing N s the scheme 
is instable. Therefore e is sufficient for estimating the convergence and stability limits. The 
measure e used here is very similar to the absolute error measure (3 for the unraveling schemes 
which has been used for studying standard first- and higher-order unraveling schemes 8 -. The 
only difference between them is the reference calculation - for calculating the error (3 one 
considers the exact solution produced numerically by some direct propagator or analytically. 
It was found 8 ' that for a small number of samples N s the error measure f3 is mainly statistical 
due to the finite sample size. For large N s the error due to the finite time step starts to 
dominate. Since both terms in the sum in Eq. (j41j) carry the same time-step error it will 
cancel out in the error measure e. Thus our convergence measure is a criterion for the 
statistical error of the stochastic scheme. 



p 1 (t j ,N s ) = ^f^Kutmm^)) +c.c] 



(40) 



As a convergence measure we introduce the quantity 




(41) 
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Considering the efficiency of the stochastic method it is convenient to look at the relation 
between numerical effort and achieved convergence. When the DDA is invoked together 
with the RWA in the system-bath coupling one can obtain the Lindblad form © that can 
be solved by the standard quantum jump metho d 4l6l8l9ll ° l17 as well as by the present one. In 
this way the numerical performance of both methods can be compared. Figure |H1 shows the 
behavior of e for k = 8 with increasing number of trajectories for the ET example with one 
reaction coordinate. Within Redfield theory the variation of T towards large system-bath 
coupling values is of limited validity. So the Redfield QME does not allow to verify the 
new scheme in the strong coupling regime. However, we performed the calculations for two 
values for T (O.luo and O.Olwo) that are supposed to belong to the weak coupling regime. In 
addition, the computation was performed for two propagation time steps St. It was found 
that for TSt « 0.1 the proposed scheme becomes unstable. As can be seen in Fig. |H] decreasing 
TSt has no influence on the convergence e. The slope of e in the double logarithmic scale 
(dloge/dlog N s ) in Fig. EI is — 1, i.e. the scheme converges as 1/N S . When simulating other, 
physically different systems we expect the convergence behavior to be not too different from 
the example studied here. At least the proportionality to 1/N a will stay unchanged. The only 
difference can be the intercept of e (i.e. the value of log epsilon for logiV s = 0) which may 
have some physical reasoning. We expect some change in the instability limit for TSt w 0.1 
with the type of system studied. This topic has to be explored in future studies, especially 
for QMEs which allow larger variation of the system-bath coupling. 

Now let us discuss the dependence on the basis size. It is well known that the numerical 
effort for solving QMEs with the quantum jump method and with direct propagators scales 
quadratically and cubically, respectively, with the basis size M . Thus, for both high accuracy 
and lower numerical expense stochastic methods would only be preferred over direct propa- 
gators if the number of trajectories N B is much smaller than M. We shall examine the scaling 
behavior by solving the one-dimensional ET model with the use of the present stochastic 
scheme for both the Redfield QME (J3j) and the Lindblad QME © as well as the standard 
scheme for the Lindblad QME (JSJ). All runs were performed with increasing basis size and 
compared to a direct propagation. For this purpose we choose the short iterative Arnoldi 
propagato r 11 ^ 2 ! 50 in a Krylov space of dimension 12. As expected, the numerical expense for 
a few trajectories is much smaller than for the direct propagator. But for converged results 
one needs much more that one trajectory. When the number of trajectories necessary for 
the complete convergence is greater than M one has to make a trade-off between the accu- 
racy achieved with the stochastic method and the numerical effort. For the comparison we 
choose as an example N s = 500 for which the convergence is not yet complete. In Fig. El 
we have seen that this sample size already yields qualitatively the same result as the direct 
propagator. For very accurate calculations one needs a much larger number of trajectories. 
The crossing point in Fig. H2 shows that for M > 212 the new stochastic algorithm should 
be preferred. In addition one gets a benefit from the efficient parallel implementation of the 
stochastic algorithm which is of great practical use especially when one is asking for fast pre- 
liminary results for the dynamics of large systems. The slope of the lines in Fig. El reflecting 
the scaling of the numerical effort for the stochastic methods is approximately 2.3 and for 
the direct propagation 3.2. These scalings are larger than the theoretical estimates of 2 and 
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3, respectively, because the true complexity of the algorithm is not a simple power function 
but a polynomial function of TV due to array operations of lower order. This difference has 
to disappear for very large Af. 

In the following the intercept of the curves for the numerical effort in Fig. El i.e. the CPU- 
time extrapolated to N = 1, is examined. The ratio of the intercepts between the proposed 
stochastic scheme for the Redfield QME and the standard scheme for the Lindblad QME is 
found to be 2.4 (see Fig.^J). The respective ratio between the proposed stochastic scheme for 
the Redfield QME and the same scheme for the Lindblad QME is 4.4. This can be interpreted 
as follows. Neglecting the operations of lower order the intercept in a double logarithmic 
plot like in Fig. El has to be proportional to the number of matrix- vector multiplications for 
one time step, i.e. the intercept holds to some extent information specific for each stochastic 
scheme. Considering only the deterministic part of the implemented algorithm (see the 
Appendix for details) each time step requires AM + 9 matrix-vector multiplications (for 
M dissipative channels, 4 x M operations for computing the jump rates Eq. (|20|l. eight 
operations for time propagation and one for calculating the population) versus only 5 for 
the standard jump algorithm (four operations for time propagation and one for calculating 
the population). Nevertheless, the values found from the simulation (Fig. EJ), 12/5 and 22/5, 
deviate from the estimates 13/5 and 17/5, respectively, due to systematic effects like lower- 
order operations for not very large M as well as due to stochastic effects like the varying 
number of quantum jumps. To summarize, for the Lindblad QME the standard algorithm 
is the method of choice, while for the Redfield QME, where one cannot apply the standard 
method, one has to cope with the larger numerical expense. The latter can be significantly 
reduced when high accuracy is not necessary. 

VII. CONCLUSION 

In this paper the stochastic unraveling technique was extended for non-Lindblad QMEs. 
This progress became possible with the use of the wave-function pair in the doubled Hilbert 
space and the derivation of stable, almost normalized SSEs. An efficient first-order quantum 
jump algorithm was proposed. The efficiency is determined by the behavior of the norm of 
every single trajectory. In this sense the jump rates were used as parameters to influence 
the efficiency. 

Occurrence of negative population for single trajectories is by no means a problem of the 
proposed unraveling scheme. Rather it is related to the fact that the QMEs (0 EJ) do not 
preserve the RDM positive semi-definite. It is known^ that the negative eigenvalues of the 
RDM in the Redfield theory arise from the inconsistency between the initial RDM and the 
bath state, i.e. due to neglected initial correlations in the Born-Markov approximation. A 
satisfactory resolution of this problem is the slippage of the initial conditions as derived by 
Gaspard et a/.—. In this method the so called slippage superoperator takes into account 
the short-time bath correlations. Applied to the initial RDM it introduces the necessary 
correlations into the initial state. This manipulation of the initial state ensures propagation 
of a positive semidefinite RDM at any further moment of time. 

The method proposed was successfully tested for the Redfield QME for the damped 
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harmonic oscillator and two ET models. It was shown that the scheme allows for more 
efficient quantum dynamical simulations of large systems. The most important benefit of 
the method is that it can be applied in a straightforward manner since the SSEs and the 
expressions of the jump rates are formally the same for other models of system and bath. 
Therefore, a potential use of the proposed method can be made in simulations involving 
any kind of non-Markovian QMEs provided that they are in a time-local form like in the 
time-convolutionless formalism 25 as well as in methods using auxiliary density matrices to 
include the memory effects^. 

Each individual trajectory occupies nearly the whole volume in phase space accessible 
for the system. It would be of possible interest to see whether the phase-space volume of 
a single trajectory shrinks with decreasing U and to study the classical behavior of the ET 
system. 
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APPENDIX: THE NUMERICAL ALGORITHM 

This algorithm gives the numerical solution of the SSEs (jl5aj) and (j!5b|) . The initial wave 
functions |^j )S (0)) and |</>j, s (0)) are constructed so that 

p s (o) = J2 w *wY, [IM°)>GM°)I + hM°)XM°)l] • ( 42 ) 

1=1 S g=l 

Here N e is the number of non-zero eigenvalues Wi of the initial density matrix p s (0) and 
N s the number of trajectories corresponding to each eigenvalue W{. The wave functions are 
propagated jointly (as pairs) as follows starting with t = 0. 

1. store/send \ipi lS (t)) and \<j>i, a {t)) for averaging; 

2. calculate the rates p\ and p\ according to Eqs. ()20a|) and (|20b[) . 

3. generate a random number e G (0, 1); 

4. ife>dt£ fc (Pfc+P!) then 
propagate \ipi, a (t)) and | 

* find \ipi 7 s{t + dt)) and \<p iiS (t + dt)) solving 

d\*Pi,s(t))/dt = A\ip iiS (t)) and d\(j) ijS (t)) / dt = A\<j> i>a (t)) , respectively 

* set t — t + dt 

* go to step [T] 
else 

* if e < dt J2k Pk tnen 
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- jump with probability p\dt: 
\il>i, s (t)) -> E k \^ iiS (t))/^/pl and \</>i,s(t)) 

else 

- jump with probability p\dt: 
\ipi,s(t)) -> C fc |V i)S (t)}/v5l and 

* go to step El 

The deterministic propagation of the wave functions in this work was performed with the use 
of a forth-order Runge-Kutta method^ 2 -. Accordingly, one time step requires four matrix- 
vector multiplication for each wave function. The ensemble-averaged expectation value of 
an observable A is calculated as 

N e N s 

P^)=E^F^ [( ^ s(t)|A| ^ ( * ))+c - c - ] • (43) 

i=l s s =l 

The method can be parallelized using MPI^. In such an implementation every single 
stochastic trajectory is propagated by a different process. Only the averaging operation 
(|4"3J) is done at certain times by means of collective communications. In this way the task 
can be efficiently distributed on a cluster of PCs. 



c*IM*)>/v®E 
Ek\Mt))/Vit 
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FIG. 1: Population dynamics of a damped harmonic oscillator simulated using the proposed 
quantum jump method with 10 3 trajectories (dashed lines) and with 10 4 trajectories (dot-dashed 
lines) compared to the direct RDM propagation (solid lines). 




FIG. 2: Relaxation of the donor population for the electron transfer model with a single reaction 
mode. The solid line shows the exact solution of the QME, the dashed line one arbitrary trajectory 
of the quantum jump method, the dotted line an average over 500 trajectories. 
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FIG. 3: Distribution of the expectation values of the population of the donor state Pi produced 
by the new unraveling scheme for the Redfield QME (dotted line) and the standard normalized 
jump method for the Lindblad QME (solid line) at time ut/(2ir) = 5.8. Both distributions are 
normalized to unity. 
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FIG. 4: Contour plots of the Wigner representation of the RDM of one individual trajectory (a), 
of the RDM recovered with 500 trajectories (b) and of the exact RDM (c). Regions with maximum 
are denoted with +, and regions with minimum with — . 
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FIG. 5: Three-dimensional plots of the Wigner representation of the RDM of one individual 
trajectory (a), of the RDM recovered with 500 trajectories (b) and of the exact RDM (c). The 
data shown here are the same as in Fig. 
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FIG. 6: Dynamics of ET depicted as a path in phase space calculated by the exact solution of 
the Redfield QME (solid line) and by averaging over 500 trajectories of the quantum jump method 
(dotted line). 




FIG. 7: Population dynamics of the excited state in the effective one-mode (thick lines) and 
the two-mode (thin lines) models for ET. Coherent dynamics are denoted by long dashed lines, 
dynamics with dissipation in DDA by dashed lines, and the Redfield dynamics by solid lines. The 
two-mode model is solved with the new stochastic method with 5000 trajectories. 



25 



10 



10 ' 



10 • 



10 



10 • 



i 9i 

*8 



'i, 



10' 



10 



10' 



10* 

N. 



10 



10" 



FIG. 8: Convergence behavior of the proposed stochastic unraveling scheme for T = O.Olwo 
(opaque), T = O.Icjq (filled), St = 1 (circles), and 5t = 10 (triangles). 
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FIG. 9: Numerical effort of the new unraveling scheme with 500 trajectories for Eq. (j3J) (circles), the 
standard quantum jump method for Eq. © (rhombs), and the new scheme for Eq. @ (triangles) 
shown for the model of one-dimensional ET. The short iterative Arnoldi method (filled squares) is 
shown as reference solving Eq. ©. The data points at M = 512 are extrapolated. 
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